home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / STATCORR.ZIP / STPRG.FOR < prev   
Text File  |  1985-11-29  |  11KB  |  324 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE STPRG
  5. C
  6. C        PURPOSE
  7. C           TO PERFORM A STEPWISE MULTIPLE REGRESSION ANALYSIS FOR A
  8. C           DEPENDENT VARIABLE AND A SET OF INDEPENDENT VARIABLES.  AT
  9. C           EACH STEP, THE VARIABLE ENTERED INTO THE REGRESSION EQUATION
  10. C           IS THAT WHICH EXPLAINS THE GREATEST AMOUNT OF VARIANCE
  11. C           BETWEEN IT AND THE DEPENDENT VARIABLE (I.E. THE VARIABLE
  12. C           WITH THE HIGHEST PARTIAL CORRELATION WITH THE DEPENDENT
  13. C           VARIABLE).  ANY VARIABLE CAN BE DESIGNATED AS THE DEPENDENT
  14. C           VARIABLE.  ANY INDEPENDENT VARIABLE CAN BE FORCED INTO OR
  15. C           DELETED FROM THE REGRESSION EQUATION, IRRESPECTIVE OF ITS
  16. C           CONTRIBUTION TO THE EQUATION.
  17. C
  18. C        USAGE
  19. C           CALL STPRG (M,N,D,XBAR,IDX,PCT,NSTEP,ANS,L,B,S,T,LL,IER)
  20. C
  21. C        DESCRIPTION OF PARAMETERS
  22. C           M    - TOTAL NUMBER OF VARIABLES IN DATA MATRIX
  23. C           N    - NUMBER OF OBSERVATIONS
  24. C           D    - INPUT MATRIX (M X M) OF SUMS OF CROSS-PRODUCTS OF
  25. C                  DEVIATIONS FROM MEAN.  THIS MATRIX WILL BE DESTROYED.
  26. C           XBAR - INPUT VECTOR OF LENGTH M OF MEANS
  27. C           IDX  - INPUT VECTOR OF LENGTH M HAVING ONE OF THE FOLLOWING
  28. C                  CODES FOR EACH VARIABLE.
  29. C                    0 - INDEPENDENT VARIABLE AVAILABLE FOR SELECTION
  30. C                    1 - INDEPENDENT VARIABLE TO BE FORCED INTO THE
  31. C                        REGRESSION EQUATION
  32. C                    2 - VARIABLE NOT TO BE CONSIDERED IN THE EQUATION
  33. C                    3 - DEPENDENT VARIABLE
  34. C                  THIS VECTOR WILL BE DESTROYED
  35. C           PCT  - A CONSTANT VALUE INDICATING THE PROPORTION OF THE
  36. C                  TOTAL VARIANCE TO BE EXPLAINED BY ANY INDEPENDENT
  37. C                  VARIABLE.  THOSE INDEPENDENT VARIABLES WHICH FALL
  38. C                  BELOW THIS PROPORTION WILL NOT ENTER THE REGRESSION
  39. C                  EQUATION.  TO ENSURE THAT ALL VARIABLES ENTER THE
  40. C                  EQUATION, SET PCT = 0.0.
  41. C           NSTEP- OUTPUT VECTOR OF LENGTH 5 CONTAINING THE FOLLOWING
  42. C                  INFORMATION
  43. C                     NSTEP(1)- THE NUMBER OF THE DEPENDENT VARIABLE
  44. C                     NSTEP(2)- NUMBER OF VARIABLES FORCED INTO THE
  45. C                               REGRESSION EQUATION
  46. C                     NSTEP(3)- NUMBER OF VARIABLE DELETED FROM THE
  47. C                               EQUATION
  48. C                     NSTEP(4)- THE NUMBER OF THE LAST STEP
  49. C                     NSTEP(5)- THE NUMBER OF THE LAST VARIABLE ENTERED
  50. C           ANS  - OUTPUT VECTOR OF LENGTH 11 CONTAINING THE FOLLOWING
  51. C                  INFORMATION FOR THE LAST STEP
  52. C                     ANS(1)- SUM OF SQUARES REDUCED BY THIS STEP
  53. C                     ANS(2)- PROPORTION OF TOTAL SUM OF SQUARES REDUCED
  54. C                     ANS(3)- CUMULATIVE SUM OF SQUARES REDUCED UP TO
  55. C                             THIS STEP
  56. C                     ANS(4)- CUMULATIVE PROPORTION OF TOTAL SUM OF
  57. C                             SQUARES REDUCED
  58. C                     ANS(5)- SUM OF SQUARES OF THE DEPENDENT VARIABLE
  59. C                     ANS(6)- MULTIPLE CORRELATION COEFFICIENT
  60. C                     ANS(7)- F RATIO FOR SUM OF SQUARES DUE TO
  61. C                             REGRESSION
  62. C                     ANS(8)- STANDARD ERROR OF THE ESTIMATE (RESIDUAL
  63. C                             MEAN SQUARE)
  64. C                     ANS(9)- INTERCEPT CONSTANT
  65. C                     ANS(10)-MULTIPLE CORRELATION COEFFICIENT ADJUSTED
  66. C                             FOR DEGREES OF FREEDOM.
  67. C                     ANS(11)-STANDARD ERROR OF THE ESTIMATE ADJUSTED
  68. C                             FOR DEGREES OF FREEDOM.
  69. C           L    - OUTPUT VECTOR OF LENGTH K, WHERE K IS THE NUMBER OF
  70. C                  INDEPENDENT VARIABLES IN THE REGRESSION EQUATION.
  71. C                  THIS VECTOR CONTAINS THE NUMBERS OF THE INDEPENDENT
  72. C                  VARIABLES IN THE EQUATION.
  73. C           B    - OUTPUT VECTOR OF LENGTH K, CONTAINING THE PARTIAL
  74. C                  REGRESSION COEFFICIENTS CORRESPONDING TO THE
  75. C                  VARIABLES IN VECTOR L.
  76. C           S    - OUTPUT VECTOR OF LENGTH K, CONTAINING THE STANDARD
  77. C                  ERRORS OF THE PARTIAL REGRESSION COEFFICIENTS,
  78. C                  CORRESPONDING TO THE VARIABLES IN VECTOR L.
  79. C           T    - OUTPUT VECTOR OF LENGTH K, CONTAINING THE COMPUTED
  80. C                  T-VALUES CORRESPONDING TO THE VARIABLES IN VECTOR L.
  81. C           LL   - WORKING VECTOR OF LENGTH M
  82. C           IER  - 0, IF THERE IS NO ERROR.
  83. C                  1, IF RESIDUAL SUM OF SQUARES IS NEGATIVE OR IF THE
  84. C                  PIVOTAL ELEMENT IN THE STEPWISE INVERSION PROCESS IS
  85. C                  ZERO.  IN THIS CASE, THE VARIABLE WHICH CAUSES THIS
  86. C                  ERROR IS NOT ENTERED IN THE REGRESSION, THE RESULT
  87. C                  PRIOR TO THIS STEP IS RETAINED, AND THE CURRENT
  88. C                  SELECTION IS TERMINATED.
  89. C
  90. C        REMARKS
  91. C           THE NUMBER OF DATA POINTS MUST BE AT LEAST GREATER THAN THE
  92. C           NUMBER OF INDEPENDENT VARIABLES PLUS ONE.  FORCED VARIABLES
  93. C           ARE ENTERED INTO THE REGRESSION EQUATION BEFORE ALL OTHER
  94. C           INDEPENDENT VARIABLES.  WITHIN THE SET OF FORCED VARIABLES,
  95. C           THE ONE TO BE CHOSEN FIRST WILL BE THAT ONE WHICH EXPLAINS
  96. C           THE GREATEST AMOUNT OF VARIANCE.
  97. C           INSTEAD OF USING, AS A STOPPING CRITERION, A PROPORTION OF
  98. C           THE TOTAL VARIANCE, SOME OTHER CRITERION MAY BE ADDED TO
  99. C           SUBROUTINE STOUT.
  100. C
  101. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  102. C           STOUT(NSTEP,ANS,L,B,S,T,NSTOP)
  103. C           THIS SUBROUTINE MUST BE PROVIDED BY THE USER.  IT IS AN
  104. C           OUTPUT ROUTINE WHICH WILL PRINT THE RESULTS OF EACH STEP OF
  105. C           THE REGRESSION ANALYSIS.  NSTOP IS AN OPTION CODE WHICH IS
  106. C           ONE IF THE STEPWISE REGRESSION IS TO BE TERMINATED, AND IS
  107. C           ZERO IF IT IS TO CONTINUE.  THE USER MUST CONSIDER THIS IF
  108. C           SOME OTHER STOPPING CRITERION THAN VARIANCE PROPORTION IS TO
  109. C           BE USED.
  110. C
  111. C        METHOD
  112. C           THE ABBREVIATED DOOLITTLE METHOD IS USED TO (1) DECIDE VARI-
  113. C           ABLES ENTERING IN THE REGRESSION AND (2) COMPUTE REGRESSION
  114. C           COEFFICIENTS.  REFER TO C. A. BENNETT AND N. L. FRANKLIN,
  115. C           'STATISTICAL ANALYSIS IN CHEMISTRY AND THE CHEMICAL INDUS-
  116. C           TRY', JOHN WILEY AND SONS, 1954, APPENDIX 6A.
  117. C
  118. C     ..................................................................
  119. C
  120.       SUBROUTINE STPRG (M,N,D,XBAR,IDX,PCT,NSTEP,ANS,L,B,S,T,LL,IER)
  121. C
  122.       DIMENSION D(1),XBAR(1),IDX(1),NSTEP(1),ANS(1),L(1),B(1),S(1),T(1),
  123.      1LL(1)
  124. C
  125. C     ..................................................................
  126. C
  127. C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
  128. C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
  129. C        STATEMENT WHICH FOLLOWS.
  130. C
  131. C     DOUBLE PRECISION D,XBAR,ANS,B,S,T,RD,RE
  132. C
  133. C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
  134. C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
  135. C        ROUTINE.
  136. C
  137. C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
  138. C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENTS
  139. C        85,90,114,132,AND 134, MUST BE CHANGED TO DSQRT.
  140. C
  141. C     ..................................................................
  142. C
  143. C        INITIALIZATION
  144. C
  145.       IER=0
  146.       ONM=N-1
  147.       NFO=0
  148.       NSTEP(3)=0
  149.       ANS(3)=0.0
  150.       ANS(4)=0.0
  151.       NSTOP=0
  152. C
  153. C        FIND DEPENDENT VARIABLE, NUMBER OF VARIABLES TO BE FORCED TO
  154. C        ENTER IN THE REGRESSION, AND NUMBER OF VARIABLES TO BE DELETED
  155. C
  156.       DO 30 I=1,M
  157.       LL(I)=1
  158.       IF(IDX(I)) 30, 30, 10
  159.    10 IF(IDX(I)-2) 15, 20, 25
  160.    15 NFO=NFO+1
  161.       IDX(NFO)=I
  162.       GO TO 30
  163.    20 NSTEP(3)=NSTEP(3)+1
  164.       LL(I)=-1
  165.       GO TO 30
  166.    25 MY=I
  167.       NSTEP(1)=MY
  168.       LY=M*(MY-1)
  169.       LYP=LY+MY
  170.       ANS(5)=D(LYP)
  171.    30 CONTINUE
  172.       NSTEP(2)=NFO
  173. C
  174. C        FIND THE MAXIMUM NUMBER OF STEPS
  175. C
  176.       MX=M-NSTEP(3)-1
  177. C
  178. C        START SELECTION OF VARIABLES
  179. C
  180.       DO 140 NL=1,MX
  181.       RD=0
  182.       IF(NL-NFO) 35, 35, 55
  183. C
  184. C        SELECT NEXT VARIABLE TO ENTER AMONG FORCED VARIABLES
  185. C
  186.    35 DO 50 I=1,NFO
  187.       K=IDX(I)
  188.       IF(LL(K)) 50, 50, 40
  189.    40 LYP=LY+K
  190.       IP=M*(K-1)+K
  191.       RE=D(LYP)*D(LYP)/D(IP)
  192.       IF(RD-RE) 45, 50, 50
  193.    45 RD=RE
  194.       NEW=K
  195.    50 CONTINUE
  196.       GO TO 75
  197. C
  198. C        SELECT NEXT VARIABLE TO ENTER AMONG NON-FORCED VARIABLES
  199. C
  200.    55 DO 70 I=1,M
  201.       IF(I-MY) 60, 70, 60
  202.    60 IF(LL(I)) 70, 70, 62
  203.    62 LYP=LY+I
  204.       IP=M*(I-1)+I
  205.       RE=D(LYP)*D(LYP)/D(IP)
  206.       IF(RD-RE) 64, 70, 70
  207.    64 RD=RE
  208.       NEW=I
  209.    70 CONTINUE
  210. C
  211. C        TEST WHETHER THE PROPORTION OF THE SUM OF SQUARES REDUCED BY
  212. C        THE LAST VARIABLE ENTERED IS GREATER THAN OR EQUAL TO THE
  213. C        SPECIFIED PROPORTION
  214. C
  215.    75 IF(RD) 77,77,76
  216.    76 IF(ANS(5)-(ANS(3)+RD))77,77,78
  217.    77 IER=1
  218.       GO TO 150
  219.    78 RE=RD/ANS(5)
  220.       IF(RE-PCT) 150, 80, 80
  221. C
  222. C        IT IS GREATER THAN OR EQUAL
  223. C
  224.    80 LL(NEW)=0
  225.       L(NL)=NEW
  226.       ANS(1)=RD
  227.       ANS(2)=RE
  228.       ANS(3)=ANS(3)+RD
  229.       ANS(4)=ANS(4)+RE
  230.       NSTEP(4)=NL
  231.       NSTEP(5)=NEW
  232. C
  233. C        COMPUTE MULTIPLE CORRELATION, F-VALUE FOR ANALYSIS OF
  234. C        VARIANCE, AND STANDARD ERROR OF ESTIMATE
  235. C
  236.    85 ANS(6)= SQRT(ANS(4))
  237.       RD=NL
  238.       RE=ONM-RD
  239.       RE=(ANS(5)-ANS(3))/RE
  240.       ANS(7)=(ANS(3)/RD)/RE
  241.    90 ANS(8)= SQRT(RE)
  242. C
  243. C        DIVIDE BY THE PIVOTAL ELEMENT
  244. C
  245.       IP=M*(NEW-1)+NEW
  246.       RD=D(IP)
  247.       LYP=NEW-M
  248.       DO 100 J=1,M
  249.       LYP=LYP+M
  250.       IF(LL(J)) 100, 94, 97
  251.    94 IF(J-NEW) 96, 98, 96
  252.    96 IJ=M*(J-1)+J
  253.       D(IJ)=D(IJ)+D(LYP)*D(LYP)/RD
  254.    97 D(LYP)=D(LYP)/RD
  255.       GO TO 100
  256.    98 D(IP)=1.0/RD
  257.   100 CONTINUE
  258. C
  259. C        COMPUTE REGRESSION COEFFICIENTS
  260. C
  261.       LYP=LY+NEW
  262.       B(NL)=D(LYP)
  263.       IF(NL-1) 112, 112, 105
  264.   105 ID=NL-1
  265.       DO 110 J=1,ID
  266.       IJ=NL-J
  267.       KK=L(IJ)
  268.       LYP=LY+KK
  269.       B(IJ)=D(LYP)
  270.       DO 110 K=1,J
  271.       IK=NL-K+1
  272.       MK=L(IK)
  273.       LYP=M*(MK-1)+KK
  274.   110 B(IJ)=B(IJ)-D(LYP)*B(IK)
  275. C
  276. C        COMPUTE INTERCEPT
  277. C
  278.   112 ANS(9)=XBAR(MY)
  279.       DO 115 I=1,NL
  280.       KK=L(I)
  281.       ANS(9)=ANS(9)-B(I)*XBAR(KK)
  282.       IJ=M*(KK-1)+KK
  283.   114 S(I)=ANS(8)* SQRT(D(IJ))
  284.   115 T(I)=B(I)/S(I)
  285. C
  286. C        PERFORM A REDUCTION TO ELIMINATE THE LAST VARIABLE ENTERED
  287. C
  288.       IP=M*(NEW-1)
  289.       DO 130 I=1,M
  290.       IJ=I-M
  291.       IK=NEW-M
  292.       IP=IP+1
  293.       IF(LL(I)) 130, 130, 120
  294.   120 DO 126 J=1,M
  295.       IJ=IJ+M
  296.       IK=IK+M
  297.       IF(LL(J)) 126, 122, 122
  298.   122 IF(J-NEW) 124, 126, 124
  299.   124 D(IJ)=D(IJ)-D(IP)*D(IK)
  300.   126 CONTINUE
  301.       D(IP)=D(IP)/(-RD)
  302.   130 CONTINUE
  303. C
  304. C        ADJUST STANDARD ERROR OF THE ESTIMATE AND MULTIPLE CORRELATION
  305. C        COEFFICIENT
  306. C
  307.       RD=N-NSTEP(4)
  308.       RD=ONM/RD
  309.   132 ANS(10)=SQRT(1.0-(1.0-ANS(6)*ANS(6))*RD)
  310.   134 ANS(11)=ANS(8)*SQRT(RD)
  311. C
  312. C        CALL THE OUTPUT SUBROUTINE
  313.       CALL STOUT (NSTEP,ANS,L,B,S,T,NSTOP)
  314. C
  315. C        TEST WHETHER THE STEP-WISE REGRESSION WAS TERMINATED IN
  316. C        SUBROUTINE STOUT
  317. C
  318.       IF(NSTOP) 140, 140, 150
  319. C
  320.   140 CONTINUE
  321. C
  322.   150 RETURN
  323.       END
  324.